source('../env.R')
Using GitHub PAT from the git credential store.
Skipping install of 'clootl' from a github remote, the SHA1 (2ed1650b) has not changed since last install.
Use `force = TRUE` to force installation
Load community data
communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))
Rows: 2428 Columns: 7── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, ebird_species_name, seasonal, presence, origin
dbl (2): city_id, relative_abundance_proxy
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
communities
Load trait data
traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_ebird.csv'))
Rows: 332 Columns: 10── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): ebird_species_name, habitat, trophic_level, trophic_niche, primary_lifestyle
dbl (5): beak_width, hwi, mass, habitat_density, migration
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(traits)
species_urban_response = communities %>% group_by(ebird_species_name) %>% summarise(mean_relative_abundance_proxy = mean(relative_abundance_proxy)) %>% left_join(traits)
Joining with `by = join_by(ebird_species_name)`
species_urban_response$is_urban = species_urban_response$mean_relative_abundance_proxy > 0
species_urban_response$ebird_species_name = str_replace(species_urban_response$ebird_species_name,' ', '_')
species_urban_response
species_urban_response %>% filter(mean_relative_abundance_proxy > 50)
species_urban_response %>% filter(mean_relative_abundance_proxy > 40)
species_urban_response %>% filter(mean_relative_abundance_proxy > 30)
species_urban_response %>% group_by(is_urban) %>% summarise(
median_beak_width = median(beak_width),
iqr_beak_width = IQR(beak_width),
median_hwi = median(hwi),
iqr_hwi = IQR(hwi),
median_mass = median(mass),
iqr_mass = IQR(mass),
n = n()
)
ggplot(species_urban_response, aes(x = is_urban, y = beak_width)) + geom_boxplot()
wilcox.test(beak_width ~ is_urban, species_urban_response)
Wilcoxon rank sum test with continuity correction
data: beak_width by is_urban
W = 1525, p-value = 0.3893
alternative hypothesis: true location shift is not equal to 0
There is not a significant response for the gape width between urban species and non urban species.
ggplot(species_urban_response, aes(x = is_urban, y = hwi)) + geom_boxplot()
wilcox.test(hwi ~ is_urban, species_urban_response)
Wilcoxon rank sum test with continuity correction
data: hwi by is_urban
W = 1213, p-value = 0.01492
alternative hypothesis: true location shift is not equal to 0
There is a significant response between urban and non-urban species
ggplot(species_urban_response, aes(x = is_urban, y = mass)) + geom_boxplot()
wilcox.test(mass ~ is_urban, species_urban_response)
Wilcoxon rank sum test with continuity correction
data: mass by is_urban
W = 1825, p-value = 0.5168
alternative hypothesis: true location shift is not equal to 0
There is not a significant response between species.
tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
tree_pruned = ladderize(drop.tip(tree, setdiff(tree$tip.label, species_urban_response$ebird_species_name)))
ggtree(tree_pruned, layout="fan")
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
columbidae_response <- species_urban_response$mean_relative_abundance_proxy
names(columbidae_response) <- species_urban_response$ebird_species_name
Is the a phylogenetic signal for our relative abundance proxy
col_ut_phylo_signal <- phylosig(tree_pruned, columbidae_response, method="lambda", test=TRUE)
col_ut_phylo_signal
Phylogenetic signal lambda : 0.46832
logL(lambda) : -576.839
LR(lambda=0) : 7.45324
P-value (based on LR test) : 0.00633224
species_urban_response$ebird_species_name_clean = gsub("_", " ", species_urban_response$ebird_species_name)
normalize <- function(x) (x - min(x)) / (max(x) - min(x)) * 5 + 1 # Scale to range [1, 6]
species_urban_response$beak_width_norm <- normalize(species_urban_response$beak_width)
species_urban_response$hwi_norm <- normalize(species_urban_response$hwi)
species_urban_response$mass_norm <- normalize(species_urban_response$mass)
g = ggtree(tree_pruned, layout="fan") %<+% species_urban_response + abundance_proxy_scale
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
g +
geom_fruit(geom = geom_point, mapping = aes(size = beak_width_norm), fill = beak_width_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0, show.legend = FALSE) +
geom_fruit(geom = geom_point, mapping = aes(size = hwi_norm), fill = hwi_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0.13, show.legend = FALSE) +
geom_fruit(geom = geom_point, mapping = aes(size = mass_norm), fill = mass_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0.15, show.legend = FALSE) +
geom_tiplab(size=2, aes(label=ebird_species_name_clean, color = mean_relative_abundance_proxy), offset = 15) +
labs(color = "Mean relative abundance proxy") +
theme(legend.position = "bottom") + ggplot2::xlim(0, 80)
ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny.jpg'), width = 2000, height = 2100, units = 'px')
cut_off_3 = 50
cut_off_2 = 35
cut_off_1 = 1
hwi_by_mass_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(hwi, mass))
hwi_by_mass_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(hwi, mass))
hwi_by_mass_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(hwi, mass))
beak_width_by_mass_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(beak_width, mass))
beak_width_by_mass_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(beak_width, mass))
beak_width_by_mass_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(beak_width, mass))
hwi_by_beak_width_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(hwi, beak_width))
hwi_by_beak_width_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(hwi, beak_width))
hwi_by_beak_width_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(hwi, beak_width))
species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% nrow()
[1] 5
species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% nrow()
[1] 15
species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% nrow()
[1] 72
polygon_line_type = 'dashed'
polygon_linewidth = 0.4
polygon_alpha1 = 0.05
polygon_alpha2 = 0.1
polygon_alpha3 = 0.15
hwi_by_mass = ggplot(species_urban_response, aes(x = hwi, y = mass, colour = mean_relative_abundance_proxy)) +
geom_polygon(data = hwi_by_mass_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = hwi_by_mass_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = hwi_by_mass_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_point() +
abundance_proxy_scale +
labs(colour = "Mean relative abundance proxy", y = 'Mass', x = '') +
theme(legend.position = 'bottom', legend.title.position = 'top', legend.key.width = unit(12, 'mm'))
beak_width_by_mass = ggplot(species_urban_response, aes(x = beak_width, y = mass, colour = mean_relative_abundance_proxy)) +
geom_polygon(data = beak_width_by_mass_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = beak_width_by_mass_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = beak_width_by_mass_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_point() +
abundance_proxy_scale +
labs(y = '', x = 'Beak width')
hwi_by_beak_width = ggplot(species_urban_response, aes(x = hwi, y = beak_width, colour = mean_relative_abundance_proxy)) +
geom_polygon(data = hwi_by_beak_width_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = hwi_by_beak_width_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = hwi_by_beak_width_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_point() +
abundance_proxy_scale +
labs(y = 'Beak width', x = 'HWI')
legend <- ggpubr::get_legend(
# create some space to the left of the legend
hwi_by_mass
)
plot_grid(
nrow = 2, ncol = 2,
hwi_by_mass + theme(legend.position="none"),
beak_width_by_mass + theme(legend.position="none"),
hwi_by_beak_width + theme(legend.position="none"),
legend
)
ggsave(filename(FIGURES_OUTPUT_DIR, 'traits.jpg'), width = 2000, height = 1800, units = 'px')